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Abstract 

A computational study of thermally-driven convection in multilayered fluid structures will be 
performed in order to examine the effect of interactions among deformable fluid-fluid interfaces on 
the structure of time-dependent flow in these systems. Multilayered fluid structures in two model 
configurations will be considered: the differentially heated rectangular cavity with a free surface, 
and the encapsulated cylindrical liquid bridge. An extension of a numerical method developed as 
part of our recent NASA Fluid Physics grant will be used to account for finite deformations of 
fluid-fluid interfaces. 


1 Introduction 

In growing crystals from a melt containing volatile components (such as GaAs or InP), stringent 
control of the stoichiometry is crucial in order to avoid crystallographic defects and degradation 
of electronic properties of the resulting product. Two approaches have been devised to minimize 
evaporation of the volatile component from the melt during processing of these materials: pres- 
surization of the growth chamber with the volatile component’s vapor, and encapsulation of the 
melt in an amorphous liquid glass phase. The former approach often requires precise tempera- 
ture control to avoid deposition of the volatile material on surfaces of the components inside the 
growth chamber. Furthermore, the pressurized vapor is often reactive (e.g., arsenic) and poses a 
safety hazard which represents a major drawback of this approach, particularly in connection with 
space-based processing. The second approach wherein the melt is encapsulated in a low melting 
point amorphous molten glass, such as boron oxide (B 2 O 3 ) or pyrolitic boron nitride (PBN), holds 
more promise for implementation in space. The viability of this technique was originally demon- 
strated by Metz et al. [25] who encapsulated germanium, PbTe, and PbSe single crystals with 
molten B 2 O 3 . Liquid-encapsulated crystal growth (abbreviated LEC growth) has since been used 
in the Bridgman and Czochralski configurations by researchers. Integration of this approach into 
the float zone process on earth has been impractical due to severe limitations on the maximum 
diameter of a stable molten zone (including the encapsulant) that can be supported by interfacial 
tension forces in the presence of destabilizing buoyancy forces. In addition, any electromagnetic 
levitation of the molten zone provided by RF heating coils would be lost upon encapsulating the 
floating zone since the encapsulant is typically nonconductive (e.g., B 2 O 3 ). Johnson [19] used a 
liquid-encapsulated float-zone (LEFZ) technique to produce 10 mm diameter ingots of polycrys- 
talline GaAs by pulling a GaAs rod through B 2 O 3 contained in a quartz tube. He estimated the 
maximum attainable crystal diameter in that process to be about 15 mm as compared with com- 
mercially useful diameters of 50 mm or more. Under microgravity conditions, however, the molten 
zone (including the encapsulant) will be affected predominantly by interfacial tension forces, and 
the limitations imposed on the maximum crystal diameter by capillary instabilities are not as se- 
vere as those dictated by buoyancy-driven instabilities on earth. In addition, there have been a 
number of interesting attempts aimed at controling interfacial instabilities through encapsulation 
of the bridge in an annular liquid column (cf. Metz et al. [25], and Barocela and Jalilevand [ 1 ]), 
and the introduction of a coflow in the encapsulating layer (cf. Lowry and Steen [23]). Hence, the 
LEFZ technique represents a promising approach for space-based production of monocrystalline 
semiconductors with a volatile component. 
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The main distinction between the unencapsulated and LEC growth techniques is the presence 
of the additional liquid-liquid interface between the melt and the encapsulant in LEC growth. 
The addition of this new interface in a highly nonisothermal environment can be a source of new 
and complex dynamics that may critically affect crystal quality. It is well-known that nonuniform 
distributions of temperature on the interface between two fluids can cause variations in interfacial 
tension across the interface. Such interfacial tension gradients act like tangential stresses on the 
interface and drive a surface flow from regions of low interfacial tension to those of high interfacial 
tension. The resulting surface flow can penetrate into the bulk phases on the two sides of the 
interface through the action of viscosity and, in turn, induce bulk fluid motions which are called 
Marangoni or thermocapillary convection. In general, thermocapillary convection will be significant 
when buoyancy effects are negligible, as in systems involving small scales or under microgravity 
conditions. 

Two of the most widely studied model problems in the analysis of thermocapillary convection 
have been the differentially-heated rectangular cavity with one free surface, and the cylindrical 
liquid bridge consisting of a liquid drop fixed between two coaxial circular disks (cf. review by 
Kuhlmann [21]). The former represents an idealization of the open-boat crystal growth technique 
while the latter is a model for the upper or lower half of the liquid domain in the float-zone process 
(cf. Ostrach [27]). For thermocapillary convection in a single fluid layer, the structure of the 
steady, linear and nonlinear, two-dimensional and axisymmetric flow and temperature distributions 
are well understood. Though significant progress has been made in understanding the nature of 
time-dependent flows, the role of interface deformations as well as the nonlinear dynamics of time- 
dependent and three-dimensional thermocapillary flows still provide challenges for future numerical 
investigations. 

In comparison with the unencapsulated process, modeling of LEC growth has received consid- 
erably less attention. Recent numerical simulations of steady thermocapillary convection in liquid 
encapsulated Czochralski growth [37, 11, 34] have predicted multicellular flow patterns in the two 
liquid phases. There have also been a few studies of steady thermocapillary convection in double 
layers of immiscible liquids within differentially heated rectangular cavities, with the imposed tem- 
perature gradient parallel to the free surface. Villers and Platten [38] performed a one-dimensional 
analysis assuming a constant temperature gradient across the cavity, while others [32, 12, 7, 22] ob- 
tained numerical solutions of the two-dimensional problem for low to moderate aspect ratios. The 
results of these studies show that under microgravity conditions, the strength of the thermocapil- 
lary flow in the layer in contact with the solid boundary can be significantly reduced, suggesting 
that the liquid encapsulation method can be used as a means of suppressing thermocapillary con- 
vection in the melt. Since the encapsulant in LEC growth is typically very viscous (e.g., B2O3), a 
reduction in the strength of the thermocapillary convection in the melt by the encapsulant is con- 
sistent with the experimental results of Eyer and Leiste [10], who showed that a solid encapsulation 
technique can eliminate striations in silicon crystals. Most recently, Prakash and Koster [30, 31] 
extended the above analyses to include a third immiscible liquid layer. They performed both a 
one-dimensional analysis similar to that of Villers and Platten [38] and an asymptotic analysis for 
shallow rectangular cavities similar to the single layer analysis of Sen and Davis [36]. The results 
of their analysis indicate that for equal layer heights and encapsulant viscosities, the flow in the 
middle (encapsulated) layer is qualitatively similar to, but weaker than, that obtained by Sen and 
Davis for a single layer. However, when the middle layer is much thicker than the encapsulant, the 
interface deformations are quite different from those for the single layer. Hence, the flow pattern 
in the encapsulated layer strongly depends on the encapsulation thickness. 

All of the above studies of thermocapillary convection in multilayered fluid structures have 
been limited to steady convection. Furthermore, aside from a few exceptions [37, 12, 31], previous 
studies have neglected deformations of the fluid-fluid interfaces in order to avoid the complexity 
associated with the determination of the unknown interface shapes. Recent numerical simulations of 
thermocapillary flows accounting for finite deformations of the free surface indicate that the shape 
of the interface may strongly influence the flow and temperature fields [24, 3, 4, 6, 20, 5, 2]. Some 
studies [28, 3, 20] have even suggested that oscillatory flows may be established as a consequence 
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of a coupling between interface deformation and thermocapillary convection. Another example of 
the importance of interface deformability in determining the dynamics of thermocapillary flows in 
LEC growth is the recent flight experiment of Koster [35] involving a three-layer fluid structure 
in a container differentially heated parallel to the free surface. The goal of the experiment was 
to identify the criteria for reducing thermocapillary convection in the middle liquid layer which is 
presumed to be representative of the encapsulated melt of an electronic material in LEC growth. 
Upon initial operation of the experiment, when two curtains separating the three fluid layers were 
drawn, the layers collapsed. Although the cause is not fully understood, the growth of capillary 
instabilities present at the fluid-fluid interfaces —as well as the contact-line region— can provide 
a plausible explanation for the failure of the experiment. 

Based on the limited available results, it is reasonable to expect the presence of multiple de- 
formable fluid-fluid interfaces in LEC growth to lead to interesting new dynamics which can play 
an important role in determining the quality of the resulting crystals grown under microgravity 
conditions. Thus, the effect of interactions among deformable fluid-fluid interfaces on the structure 
of time-dependent thermocapillary convection in these systems remains an open question that war- 
rants further investigation. It is the purpose of this study to address this issue for thermocapillary 
convection in multilayered fluid structures in the cavity and bridge configurations. We intend to 
remove the limitations associated with the assumption of nondeformability of interfaces in previous 
studies of this problem by using an extension of a numerical method developed as part of our re- 
cent NASA Fluid Physics grant [15, 16] to allow interface shapes to evolve in accordance with the 
physics of the problem. The results will help extend our knowledge of the dynamics and stability of 
multilayered fluid systems, and will provide insight into the governing physics of the LEC growth 
process. 

2 Formulation 

In what follows, the subscript ‘s’ denotes an interface-specific quantity, while a plain symbol or one 
with a caret references a property of the bulk phase on one or the other side of the interface. The 
interfacial tension, < 7 , is taken to depend on the temperature, T, as well as the surface concentration 
of surface-active contaminants, T. The temperature dependence of other fluid properties can also 
be easily included. Using the velocity scale U = (GL/fi)a tT , obtained through a tangential stress 
balance at an interface, the following well-known parameters arise in the equations governing the 
various transport processes in each phase: 

Re=UL/u , Ma = UL/a, Ca = nU/a, Bi = hL/k, Br = (nU 2 /kGL), 

where G denotes the imposed temperature gradient and h is the heat-transfer coefficient at the 
liquid-vapor interface. In addition to the above dynamic and thermal parameters, there are also a 
number of geometric and physical property ratios for the various fluid phases. The vector-invariant 
forms of the governing equations in each fluid phase are given by 

u, t + V(uu) = ^v n, V • u = 0, T,t + V-(uT) = [v 2 T + Brr : Vu] , 

where II = —pi + r is the Newtonian stress tensor with r = Vu + (Vu)'. 

The interfacial boundary conditions on x = x s (t) include 

n • (II - An) + V 3 <7 = a(V, • n)n, u = u, dxjdt - u, T = T , n • VT = (3n • VT , 

where n denotes the unit normal vector at the interface. The thermal boundary condition at the 
liquid-vapor interface is -n VT = Bi(T - T 0 ) for some ambient temperature T 0 . The solid-liquid 
boundary conditions are standard, except in the vicinity of the contact line, where a multi-valued 
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velocity gives rise to a stress singularity (cf. Dussan V. and Davis [8]). This singularity will be 
regularized through the introduction of a slip coefficient, a, according to (Huh and Scriven [18]) 

b ■ IT * e = au • e, (a«1), 

where b and t are the unit vectors normal to the solid-liquid boundary and tangent to the contact 
line, respectively, and e = t A b. 

In the presence of bulk-insoluble surface contamination, the concentration T of the surface- 
active impurity on a fluid-fluid interface is governed by 

^ + v s • (u 4 r) + (u • n)(v, • n)r = -^-v?r , 

where Pe s denotes a surface Peclet number. This surface convective-diffusion equation will be 
converted into a volumetric one and treated as such (cf. [17]): 

3 Method 

As part of a previous MSAD-funded effort, a computer code for the simulation of flows driven, or 
affected, by interfacial forces was developed; its main features are: 

* Solution of the three-dimensional, incompressible Navier-Stokes equation, 

* Temperature-dependent properties, 

* Automatic adaptive grid refinement near fluid-fluid interfaces, 

* Determination of surface topology (e.g. n, curvature, etc.), and normal distance to the 
interface via the solution of the level-set 1 equations (cf. Osher and Sethian [26]), 

* Smearing of the interfacial forces into a body force occupying a thin band about the original 
crisp interface, via a mollified delta function (cf. Peskin [29]) in the interface-normal direction. 

This code has been well calibrated and tested, and is being parallelized in order to achieve sig- 
nificantly faster speeds. The adaptive local refinement of the underlying grid in the vicinity of 
fluid-fluid interfaces is an essential feature of the code (cf. Haj-Hariri et al. [14]), and has been 
achieved using an oct-tree data structure. 

Briefly, the computational approach involves the treatment of interfacial forces as body forces, 
and smearing of all fluid properties across the interface, resulting in an effectively single-phase flow 
problem. For example, the momentum equation involves a body-force concentration, f, accounting 
for the effect of an interface: 


f = [(I - nn) • Vcr — ( V * n)(7n] b(<p) . 

In the above expression, 6 denotes a mollified Dirac delta function (cf. Peskin [29]) and <£ is the 
’distance’ from the interface, obtained via the level-set construction of Osher and Sethian [26]: 

|V0| = 1, <ps = 0 . 

The gaussian curvature of the interface ,V ■ n, is then obtained using n = V<j>. The treatment of the 
surface convective-diffusion equation (arising when an interface is contaminated with surface-active 
impurities) is more challenging, as the physical quantity P is genuinely a surface quantity and its 
treatment as a volume quantity generates ambiguity for its behavior in the normal direction. To 
avert the ambiguity, it is shown that T should be extended away from the interface using 

n ■ VT = 0 , 

*The concept of the level sets is not new; it is akin to the constant-phase surfaces which are obtained via the 
solution of the Eikonal equation in monochromatic wave-propagation problems 
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so that the surface convective-diffusion equation is then equivalent to a volumetric one. 

Another point requiring special treatment is the contact-line area. The slip-coefficient formula- 
tion is used to regularize the force singularities. However, one needs to further implement a model 
accounting for the contact-angle hysteresis as well as the general dependence of the advancing and 
receding contact angle on the in-plane normal speed of the contact line (e.g. power-law dependence 
jj — c f. Ehrhard and Davis [9]). This information is incorporated into the formulation 

through the imposition of the interface normal on the interface at the contact line, i.e. n is deter- 
mined such that nb = cos# (while n*t = 0). This specification supercedes the definition n = V<£, 
which is applied at all other points within the computational domain. 

While the stability problem will be addressed via a complete time-dependent solution, a linear 
stability analysis in the spirit of Ramanan and Homsy [33] will also be performed in order to 
illuminate the energetics of the flow. Namely, any physical quantity, q , will be expanded about it’s 
value at time t: 

q(t + St) = q(t) + Sq , 

and a linearized equation for Sq, as well as an evolution equation for an energy-like quadratic, will 
be determined. The latter equation will provide further insight into the nature of the otherwise 
purely-numerical solution by indicating the main contributions to the production and dissipation 
of energy; this information will be helpful in understanding and controlling the behavior of the 
solution. 


References 

[1] E. Barocela and A. Jalilevand. Liquid encapsulated float zone method for microgravity pro- 
duction of galium arsenide. Paper 87-0390, AIAA aerospace sciences meeting, Reno, 1987. 

[2] C. Chen and J. M. Floryan. On the thermocapillary convection in liquid layers. In Fluid 
Mechanics Phenomena in Microgravity, edited by D. A. Siginer, R. L. Thompson, and L. M. 
Trefethen, ASME, 1993. 

[3] J. C. Chen and F. S. Hwu. Oscillatory thermocapillary flow in a rectangular cavity. Int. J. 
Heat Mass Transfer , 36(15):3743-3749, 1993. 

[4] J. C. Chen, J. C. Sheu, and S. S. Jwu. Numerical computation of thermocapillary convection 
in a rectangular c avity. Numer. Heat Transfer A, 17:287-308, 1990. 

[5] J. C. Chen, J. C. Sheu, and Y. T. Lee. Maximum stable length of nonisothermal liquid bridges. 
Phys. Fluids A, 2(7):1118-1123, 1990. 

[6] C. Cuvelier and J. M. Driessen. Thermocapillary free boundaries in crystal growth. J. Fluid 
Mech., 169:1-26, 1986. 

[7] T. Doi and J. N. Koster. Thermocapillary convection in two immiscible liquid layers with free 
surface Phys. Fluids A, 5(8): 1914—1927, 1993. 

[8] E. B. Dussan V. and S. H. Davis. On the motion of a fluid-fluid interface along a solid surface. 
J. Fluid Mech., 65:71-95, 1974. 

[9] P. Ehrhard and S. H. Davis Non-isothermal spreading of liquid drops on horizontal plates J. 
Fluid Mech., 229:365, 1991 

[10] A. Eyer and H. Leiste. Striation-free silicon crystals by float-zoning with surface-coated melt. 
J. Crystal Growth, 71:249-252, 1985. 

[11] J. P. Fontaine, A. Randriamampianina, G. P. Extremet, and P. Bontoux. Simulation of steady 
and time-dependent rotation-driven regimes in a liquid-encapsulated Czochralski configuration. 
J. Crystal Growth, 97:116-124, 1989. 

[12] J. P. Fontaine and R. L. Sani. Thermocapillary effects in a multilayered fluid system, paper 
92-0689, AIAA aerosapce sciences meeting, Reno, NV, 1992. 

[13] H. P. Greenspan and B. M. McCoy. On the wetting of a surface by a very viscous fluid. Studies 
in Appl. Math,, 64:95-12, 1981. 

[14] H. Haj-Hariri, Q. Shi, and A. Borhan. Effect of local property smearing on global variables : 
Implication for numerical simulations. Phys. Fluids A, 6(8):2555 2557, 1994. 


293 



[15] H. Haj-Hariri, Q. Shi, and A. Borhan. Thermocapillary motion of deformable drops at finite 
Reynolds and Marangoni numbers. J. Fluid Mech., submitted. 

[16] H. Haj-Hariri, Q. Shi, and A. Borhan. A continuum surface model for flows driven by inter- 
facial tension gradients. Proceedings of the 17th Boundary Element International Conference , 
Madison, WI, 1995. 

[17] H. Haj-Hariri and A. Borhan. A computational model for surfactant transport on deformable 
two-phase boundaries based on interface mollification. To appear in Proceedings of the 18th 
Boundary Element International Conference, Braga, Portugal, 1996. 

[18] C. Huh and L. E. Scriven. Hydrodynamic model of steady movement of a solid/liquid/fluid 
contact line. J. Colloid Interface Sci ., 35:85-101, 1971. 

[19] E. S. Johnson. Liquid encapsulated float zone melting of GaAs. J. Crystal Growth , 30:249-256, 
1975. 

[20] N. D. Kazarinoff and J. S. Wilkowski. Bifurcations of numerically simulated thermocapillary 
flows in axially symmetric float zones. Phys. Fluids A , 2:1797-1807, 1990. 

[21] H.C. Kuhlmann. Thermocapillary flows in finite-size systems. Math . and Computer Modelling , 
20( 10-1 ): 145-173, 1994. 

[22] Q. S. Liu, G. Chen, and B. Roux. Thermogravitational and thermocapillary convection in 
a cavity containing two superposed immiscible liquid layers. Int. J. Heat Mass Transfer , 
36(1):101-117, 1993. 

[23] B. J. Lowry and P. H. Steen. Stabilization of an axisymmetrical liquid bridge by viscous- flow. 
Int, J. Multiphase Flow , 20(2):439-443, 1994 

[24] W. Q. Lu. Boundary element analysis of thermocapillary convection with a free surface in a 
rectangular cavity. Int. J. Heat Mass Transfer , 37(7): 1063—1071 , 1994. 

[25] E. P. Metz, R. C. Miller, and R. Mazelsky. A technique for pulling single crystals of volatile 
materials. J . Appi Phys., 33:2016-2017, 1962. 

[26] S. Osher and J. A. Sethian. Fronts propagating with curvature-dependent speed: Algorithms 
based on hamilton-jacobi formulations. J . Comp. Phys., 79:12-49, 1988. 

[27] S. Ostrach. Low-gravity fluid flows. Ann. Rev. Fluid Mech., 14:313-345, 1982. 

[28] S. Ostrach, Y. Kamotani, and C. L. Lai. Oscillatory thermocapillary flows. PCH PhysicoChem. 
Hydrodyn., 6:585-599, 1985. 

[29] C. S. Peskin. Numerical analysis of blood flow in the heart. J. Comp. Phys., 25:220-252, 1977. 

[30] A. Prakash and J. N. Koster. Natural and thermocapillary convection in three layers. Eur. J. 
Mech., B/Fluids, 12(5):635-655, 1993. 

[31] A. Prakash and J. N. Koster. Convection in multiple layers of immiscible liquids in a shallow 
cavity. 1. steady natural-convection. Int. J. Multiphase Flow , 20(2):383-396, 1994. 

[32] N. Ramachandran. Thermal buoyancy and Marangoni convection in a two fluid layered system- 
a numerical study. Paper 90-0254, AIAA aerospace sciences meeting, Reno, NV, 1990. 

[331 N. Ramanan and G. M. Hornsy. Linear stability of lid-driven cavity flow. Phus. Fluids , 
6(8):2690-2701, 1994. 

[34] P. Sabhapathy and M. E. Salcudean. Numerical analysis of heat transfer in LEC growth of 
GaAs. J. Crystal Growth, 97:125-135, 1989. 

[35] G. Seibert. Early mission report on the four ESA facilities: Biorack; Bubble, Drop and Particle 
Unit; Critical Point Facility and Advanced Protein Crystallization Facility flown on the IML-2 
Spacelab Mission. Microgravity News from ESA, 7(3):2-7, 1994. 

[36] A. S. Sen and S. H. Davis. Steady thermocapillary flows in two-dimensional slots. J . Fluid 
Mech., 121:163-186, 1982. 

[37] P. D. Thomas, J. J. Derby, L. J. Atherton, R. A. Brown, and M. J. Wargo. J. Crystal Growth, 
96:135-152, 1989. 

[38] D. Villers and J. K. Flatten. Influence of interfacial tension gradients on thermal convection 
in two superposed immiscible liquid layers. Appl. Sci. Res., 47:177-191, 1990. 


294 



